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Abstract. The pencil-beam model is valid only when elementary Gaussian 
beams are small enough with respect to lateral heterogeneity of a medium, which is 
not always the case in heavy charged particle radiotherapy. This work addresses 
a solution for this problem by applying our discovery of self-similar nature of 
Gaussian distributions. In this method, Gaussian beams split into narrower and 
deflecting daughter beams when their size has exceeded the lateral heterogeneity 
limit. They will be automatically arranged with modulated areal density for 
accurate and efficient dose calculations. The effectiveness was assessed in an 
carbon-ion beam experiment in presence of steep range compensation, where 
the splitting calculation reproduced the detour effect of imperfect compensation 
amounting up to about 10% or as large as the lateral particle disequilibrium effect. 
The efBciency was analyzed in calculations for carbon-ion and proton radiations 
with a heterogeneous phantom model, where the splitting calculations took about 
a minute and were factor of 5 slower than the non-splitting ones. The beam- 
splitting method is reasonably accurate, efficient, and general so that it can be 
potentially used in various pencil-beam algorithms. 



PACS numbers: 87.55. D-,87.55.Kd 

1. Introduction 

In treatment planning of radiotherapy with protons and heavier ions, the pencil-beam 
(PB) algorithm is commonly used (Hong et al 1996, Kanematsu et al 1998, 2006, 
Schaffner et al 1999, Kramer et al 2000), where a radiation field is approximately 
decomposed into two-dimensionally arranged Gaussian beams that receive energy loss 
and multiple scattering in matter. In the presence of heterogeneity, these beams grow 
differently to reproduce reahstic fluctuation in the superposed dose distribution. 

Comparisons with measurements and Monte Carlo (MC) simulations, however, 
revealed difficulty of the PB algorithm at places with severe lateral heterogeneity such 
as steep areas of a range compensator and lateral interfaces among air, tissue, and 
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bone in a patient body (Goitein 1978, Petti 1992, Kohno et al 2004, Ciangaru et al 
2005). One reason for the difRculty is that particles in a pencil beam are assumed 
to receive the same interactions, whereas they may be spatially overreaching beyond 
the density interface. The other reason is that only straight paths radiating from a 
point source are considered in beam transport, whereas actual particles may detour 
randomly by multiple scattering. 

Schneider et al (1998) showed that a phase-space analysis could address the 
overreach and detour effects for a simple lateral structure. Schaffner et al (1999) and 
Soukup et al (2005) subdivided a physical spot beam virtually into smaller beams to 
naturally reduce overreaches. Pflugfelder et al (2007) quantified lateral heterogeneity, 
with which subdivision and arrangement could be optimized. Unfortunately, those 
techniques are ineffective against beam-size growth during transport. 

For electrons, the overreach and detour effects are intrinsically much severer. 
Shiu and Hogstrom (1991) developed a solution, the PB-redefinition algorithm, 
where minimal pencil beams are occasionally regenerated, considering electron flows 
rigorously. The same idea was in fact partly applied to heavy particles for beam 
customization (Kanematsu et al 2008b), but the poly-energetic beam model to 
deal with heterogeneity could be seriously inefficient in high-resolution calculations 
necessary for Bragg peaks. 

In this study, we develop an alternative method to similarly address the overreach 
and detour effects. In the following sections, we incorporate our findings on the 
Gaussian distribution into the PB algorithm, test the new method in a carbon-ion 
beam experiment, and discuss the results and practicality for clinical applications. 



2. Materials and methods 



2.1. Theory 

2.1.1. Pencil-beam algorithm The PB algorithm in this study basically follows our 
former works (Kanematsu et al 1998, 2006, 2008b). A pencil beam with index b 
is described by position rf,, direction Vb, number of particles rib, residual range Rb, 
angular variance 0^f,, angular-spatial covariance 9tb, and spatial variance i^f, of the 
involved particles. As described in Appendix A these parameters are initialized and 



modified with transport distance s. The resultant beams with variance cr^ = t'^b are 
superposed to form dose distribution 

sr-^ Ub D^o{dbr) f \fQb + SbrVb-r\'^\ 

- ? [ R(^ j ' 

Sbr = (r- fob) ■ Vb, dbr = Rq - Rbisbr), (2) 

where fob is the beam-6 origin, Sbr is the distance at the closest approach to point r, 
dbr is its equivalent water depth, and I?$o and Ro are the tissue-phantom ratio and 
the beam range in water. 



2.1.2. Self- similarity of Gaussian distribution Any normalized Gaussian distribution 
Gm.a-{x) with mean m and standard deviation a can be represented with the standard 
normal distribution N{x) = Go,i{x) as 

G„.,.(x) = -ivf^^V N{x)^^e-^'/\ (3) 
V / v27r 
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Figure 1. The standard normal distribution N{x) (gray area) and its 
approximate functions (a) N2(x), (b) N3{x), and (c) Nji(x) (solid lines) comprised 
of multiple, displaced, narrowed, and scaled Gaussian distributions (dashed lines). 



Table 1. Size-reduction, displacement, and share-fraction factors for Gaussian 
function splitting of multiplicity m. 



Factor name 


Symbol 


m = 2 


m = 


3 


m = 4 


Size reduction 
Displacement 
Share fraction 


Tm 
dm 
fm 


i) 

V 2 ' 2 / 


1 

("1,0, 
(- - 

V 4 ' 2 ' 


+1) 
\) 


1 

2 

I -3 -1 +1 -l-3\ 
\ 2 ' 2 ' 2 ' 2 / 
(•13 3 IN 
V8 ' 8 ' 8 ' 8/ 



Incidentally, wc have found that binomial Gaussian funetion 



N2{x) = \ 



G 



n/3 



(4) 



reasonably approximates N(x) as shown in figure [TJa), where we first fixed symmetric 
displacement to = ±1/2 for the binomial terms and determined their reduced standard 
deviation a = \/3/2 to conserve variance N2{x) da; = 1. Similarly, the daughter 

Gaussian terms in N2{x) splits into grand daughters to form approximate function 



(5) 



and then into grand-grand daughters to form approximate function 

G'_3.i(x) -l-3G_i +3Gi i(x) -l-G| 



Ni{x) = I ^ 



as shown in figures [Ijb) and[ljc). Table [1] summarizes size-reduction, displacement, 
and share-fraction factors for splitting with Nm (to e {2, 3, 4}). Further splitting with 
the same displacement is not possible with valid (cr > 0) Gaussian terms. 

An overreaching Gaussian beam may split two-dimensionally into to x to smaller 
beams with these approximations. Because beam multiplication will explosively 
increase computational amount, it must be applied only when and where necessary 
with optimum multiplicity to for required size reduction. 



2.1.3. Lateral heterogeneity In a grid- voxel patient model with density distribution 
ps(f), we define density gradient vector Vps as 

VPS = 2^ 7 eg, (7) 
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Figure 2. Definitions of the source coordinate system {x,y,z), the beam-6 
coordinate system {s,t,u), and the density/dose grids at spacing Si, 82, and <53 
along axes ei, 62, and 63. The axes normal to the viewing plane, eu and 63, are 
not shown. 



where Sg and eg are the grid interval and the basis vector for axis g G {1,2,3} as 
shown in figure [2] and operation maxa[a, 6] equals a if \a\ > \b\ or otherwise b. We 
quantify the lateral heterogeneity by effective lateral density gradient 



with which we define the distance to an interface of density change Kp as 

dint(f) = minf^^, 2 S^y] , (9) 

\lxy[r) / 

where Kp — 0.1 may be appropriate for interfaces among air {ps ~ 0), soft tissues 
(0-9 < PS < 1-1), and bones (1.2 < ps < 1-7) (Kanematsu et al 2003). Effective 
lateral grid interval 

3 

<5vox = ^<5g Cg, (10) 

multiplied by 2 is the effective distance to a second laterally adjacent grid, beyond 
which the distance to the interface can not be estimated from the gradient. 

2.1.4- Beam splitting The pencil beams are examined at every transport step in a 
patient. Ones subject to splitting should be not only overreaching beyond a density 
interface but also substantially influential to the dose for computational efficiency. We 
thus require conditions 

nb>Knnbo, Rb>KRRo, (7b>-%, o-fc > dint, (11) 

V6 
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for mother beam b to spht. Cutoff parameters k„ = 0.1 and kji = 0.1 set limits on 
number of particles and residual range Ri, with respect to those of the ancestral 
original beam, rifto and Rq- Condition ab > Sr^y/y/E suppresses splitting into beams 
narrower than effective grid resolution Sxy/V^ with m — 2. Condition ai, > dint 
defines the state of overreaching. The optimum multiplicity is determined as 

2 for (T2 (Tb < dint 

3 for (73 O-f, < dint < 0-2 CTb (12) 

4 for (73 o-f, > dint, 

to suppress recursive splitting with m — 2 and 3. With the beam-6 coordinate system 
(s, i, u) shown in figure [2] and defined as 

s ^ {r- fob) ■ es, t ^ {f- fob) ■ et, u = {f - f^b) ■ e^, (13) 
es^Vb, et = — , „ „ , e„ = X et, (14) 

daughter beams b'^^ {a, (3 G [1, m]) are initialized as 

^Ki3 =^b + <yb {d„ip et + dma e«) , (15) 



fb'^^ -^b + =^fb I fb'^„ ^^b + =^Vb] , (16) 




nb'^^ — fma fmpnb, Rb'^^^ — Rb, t'^b'^^^ — (^m.t'^b, (17) 



9tb'^^=aldtb, 9\^^dh-{l-<7l)^, (18) 

i b 

where fb' ^ is the displaced position, vy ^ is the radial direction from the focus 
or the virtual source (ICRU-35 1984) of the mother beam, ny ^ is the number of 

shared particles, Rb'^^ is the conserved residual range, f^y^^ is the reduced spatial 
variance, and Oty ^ and O'^y ^ conserve focal distance {t'^/6t) and local angular variance 

(02 — Ot' /t'^) in splitting. The mother beam splits into the daughter beams to form 
different detouring paths. 

Sets of the initial parameters for daughter beams are sequentially pushed on the 
stack of computer memory and the last set on the stack will be the first beam to 
be transported in the same manner, which will be repeated until the stack has been 
emptied before moving on to the next original beam. 



2.2. Experiment 

2.2.1. Apparatus An experiment to assess the present method was carried out with 
accelerator facility HIMAC at National Institute of Radiological Sciences. A i2q6+ 
beam with nucleon kinetic energy E/A = 290 MeV was broadened to a uniform field 
of nominal 10-cm diameter by the spiral- wobbling method (Yonai et al 2008). The 
horizontal wobbler at z = zx = 527 cm and the vertical wobbler at z = zy = 470 cm 
formed a spiral orbit of maximum 10-cm radius on the isocenter plane. A 0.8-mm-thick 
Pb (ps = 5.77, Xo — 0.561 cm) foil was placed at z = 425 cm as a scatterer, which 
increased the instantaneous RMS beam size from pristine 8.3 mm to 25 mm at the 
isocenter. A large-diameter parallel-plate ionization chamber was placed at z ~ 400 
cm for dose monitoring and beam-extraction control. An Al (ps = 2.12, Xq = 8.90 
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Figure 3. Experimental layout of the beam delivery, customization, and phantom 
systems. 



cm) ridge filter for semi-Gaussian range modulation of m = 0.54 cm and a — 0.18 
cm in water (Schaffner et al 2000) and a 2-mm-tliick Al base plate were inserted at 
z = 235 cm to moderate the Bragg peak just to ease dosimetry. 

As shown in figure [31 a water (ps = 1, -'^o = 36.08 cm) tank with a 1.9-cni-thick 
PMMA (ps — 1.16, Xq — 34.07 cm) beam-entrance wall was placed at the irradiation 
site with the upstream face at z = 16.9 cm. The radiation field was defined by a 
8-cm-square 5-cm-thick brass collimator whose downstream face was at 65 cm. Two 
identical 3-cni-thick PMMA plates were inserted. The downstream plate was attached 
to the beam-entrance face of the tank covering only the x > Q side to form a phantom 
system with a bump. The upstream plate was put with its downstream face at z = 35 
cm covering only the a; < side to compensate the bump. Such arrangement is typical 
for range compensation and sensitive to the detour effects (Kohno et al 2004). These 
beam-customization elements were manually aligned to the nominal central axis at an 
uncertainty of 1 mm. 

A multichannel ionization chamber (MCIC) with 96 vented sense volumes aligned 
at intervals of 2 mm along the x axis was installed at y = in the water tank. The 
MCIC system was electromechanically movable along the z axis and the upstream limit 
at Zicf = 14.77 cm was chosen for the reference point with reference depth dj-cf = 2.44 
cm of equivalent water from the tank surface. 



2.2.2. Measurement With a reference open field without the PMMA plates or the 
collimator, we measured reference dose/MU reading Mrefi at reference height Zrof for 
every channel i for a calibration purpose. Every dose/MU reading Mi{z) of channel i 
at height z for any field is divided by corresponding reference reading M-^cH to measure 
dose D at position [xi , z) as 

Dix..z)^^^ (19) 

Airofi Zx — Zrcf ZY — Z^cf 

where divergence-correction factor zx/{zx — ^^ref ) • zy / {zy — Zrcf ) = 1-054 is to measure 
the doses in dose unit U that would be the isocenter dose for the reference depth of 
the reference field. 

We then measured reference-field doses Do{z) in the phantom at varied z 
positions, from which we get tissue-phantom ratio 

Zx — Z Zy — Z 

D^o(d) = Dq[Z) , d = drcf + Zrof - 2- (20) 

ZX Zy 
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Figure 4. (a) Tissue-phantom ratio D^o{d) witli indications for the measurement 
(o) and the reference and 80%-dose depths (dref and dgo) and (b) effective density 
ps{x,z) and (c) effective lateral density gradient ■yj:y{x,z) distributions in gray 
scale at y = in the calculation model. 



Table 2. Estimated contributions of beam-line elements at height z to beam 
range (R), scattering angle (8), and source sizes (cx and cry). 















Element 


z 


-AR 








Pristine 




0.68 cm 




8.3 mm 


8.3 mm 


Scatterer 


425 cm 


0.46 cm 


5.5 mrad 


5.6 mm 


2.5 mm 


Ridge filter 


235 cm 


0.96 cm 


3.2 mrad 


9.3 mm 


7.5 mm 


Total 




2.10 cm 




13.7 mm 


11.5 mm 



Beam range Rq with Gaussian modulation was equated to the distal 80%-dose depth 
(Koehler et al 1975) dso = 14.14 cm as shown in figure [D^a) . 

With the collimator and the PMMA plates in place, lateral dose profiles were 
measured in the same manner with particular interest around z = 3.3 cm, 6.8 cm, and 
10.3 cm, where the Bragg peaks were expected for the primary ions passing through 
none, either, and both of the PMMA plates. 

2.2.3. Calculation Table [2] shows range loss ~AR and scattering A9^ for the beam- 
line elements, and the resultant contributions to source sizes ax and ay estimated by 
back projection to the sources. The ridge filter with the base plate was modeled as 
plain aluminum of average thickness. The scattering for the scatterer was estimated 
from measured beam size 25 mm quadratically subtracted by pristine size in the 
distance of 425 cm. Total range loss 2.10 cm was deduced from range 16.24 cm 
expected for E/A = 290 MeV carbon ions (Kanematsu 2008c) and deficit 0.68 cm for 
the pristine beam may be attributed to minor materials in the beam line. 

As described in [Appendix A[ pencil beams were defined to cover the coUimated 
field at intervals of 6i — I mm on the isocenter plane, where the open field was assumed 
to have uniform unit fluence $o = ni,/S'^ — 1. Exact collimator modeling was omitted 
because we were interested in the density interface in the middle of the field. The 
upstream PMMA plate was modeled as a range compensator with range loss S = 3.48 
cm for a; < or S' = for x > 0, where the original beams were generated, followed 
by the range loss and scattering. The phantom system comprised of the downstream 
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PMMA plate and the water tank was modeled as density voxels at grid intervals of 
Si = S2 = 5^ = 1 mm for a 2-L volume of < 5 cm, \y\ < 5 cm, and < z < 20 cm. 
Figures Ufb) anddjc) show the density and lateral heterogeneity distributions. 

We carried out dose calculations with beam splitting enabled (splitting 
calculation) and disabled (non-splitting calculation). In this geometry, the density 
interface at a; = was almost parallel to the beams and only ones in the two nearest 
columns would split. 

2.3. Applications 

To examine effectiveness and efficiency of this method with larger heterogeneity, a 3- 
cm diameter cylindrical air cavity at (x, z) = (—3 cm, 13 cm) and two 1-cm diameter 
bone rods with density ps = 2 at (2 cm, 13 cm) and (4 cm, 13 cm) were added to the 
phantom in the calculation model. 

We carried out splitting and non-splitting dose calculations of the same carbon- 
ion radiation to monitor changes in frequencies of splitting modes, number of 
stopped beams, total path length J2b I total effective volume J2b 1 12(T^ds in 
the heterogeneous phantom, and computational time with a 2-GHz PowerPC G5/970 
processor by Apple/IBM. 

The splitting calculation could be more effective for protons because they 
generally suffer larger scattering. We thus carried out equivalent dose calculations 
for protons with enhanced scattering angle by factor 3.61 (jA.6|) in otherwise the same 
configuration including the tissue-phantom-ratio data. 

3. Results 

3.1. Experiment 

Figure [5] shows the two-dimensional dose distributions measured in the carbon- 
ion beam experiment and the corresponding non-splitting and splitting calculations. 
Figure [6] shows their lateral profiles in the plateau and at depths for sub peak, main 
peak, and potential sub peak expected for particles that penetrated both, either, and 
none of the PMMA plates. A dip/bump structure was commonly formed along the 
X = line for lateral particle disequilibrium (Goitein 1978). There was actually a 
sub peak in the measurement and in the splitting calculation, while it was naturally 
absent in the non-splitting calculation. The observed loss of the main-peak component 
was also reproduced by the splitting calculation. The potential sub peak was barely 
noticeable only in the splitting calculation. 

3.2. Applications 

Figure [7] shows details of the heterogeneous phantom and the dose distributions by 
splitting calculation for the carbon-ion and proton radiations. The larger scattering 
for protons naturally led to the larger dose blurring. Figure [5] shows the dose profiles 
at the main peak and where the heterogeneity effects were large by splitting and 
non-splitting calculations. In addition to the loss of the main-peak component at 
a; « 0, beam splitting caused some dose enhancement in the shoulders of the profiles 
especially for the carbon ions. 
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Figure 5. Dose distributions in linear gray scale at 2/ = by (a) measurement, 
(b) splitting calculation, and (c) non-splitting calculation. 
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Figure 6. Lateral dose profiles at y = by measurement (o), splitting calculation 
(solid), and non-splitting calculation (dashed) at (a) z = 14.8 cm (plateau), (b) 
10.3 cm (sub peak), (c) 6.8 cm (main peak), and (d) 3.3 cm (potential sub peak). 
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Figure 7. Distributions in the heterogeneous phantom at y = of (a) density pg 
and (b) effective lateral density gradient 'yxy in the calculation model and doses 
D/U from (c) carbon-ion and (d) proton radiations calculated with splitting. 
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Figure 8. Lateral dose profiles in the heterogeneous phantom by splitting 
(solid) and non-splitting (dashed) calculations for projectile/height of (a) carbon 
ion/7.3 cm, (b) carbon ion/6.8 cm (main peak), (c) carbon ion/5.8 cm, (d) 
proton/7.3 cm, (e) proton/6.8 cm (main peak), and (f) proton/5.8 cm. 



Table 3. Statistics per original beam in splitting and non-splitting calculations 
for carbon-ion and proton beams in the heterogeneous phantom. 



Projectile 


Carbon ion 


Proton 


Beam splitting 


No 


Yes 


No 


Yes 


Frequency of m = 2 





0.243 





3.813 


Frequency of m = 3 





0.132 





0.714 


Frequency of m = 4 





1.636 





0.967 


Number of stopped beams 


1 


26.8 


1 


25.0 


Mean path length/cm 


20.0 


394.6 


20.0 


499.8 


Mean effective volume/cm'^ 


3.52 


23.1 


30.8 


380.4 


Computational time/s 


9.3 


45.8 


15.3 


63.8 



Table [3] shows the statistical results, where the splitting effectively increased the 
carbon-ion and proton beams by factors of 27 and 25 in number, 20 and 25 in path 
length, 6.6 and 12 in volume, and 4.9 and 4.2 in total computation. 

4. Discussion 

Subdivision of a radiation field into virtual pencil beams is an arbitrary process in 
the PB algorithm although the beam sizes and intervals should be limited by lateral 
heterogeneity of a given system. In the PB-redefinition algorithm (Siu and Hogstrom 
1991), beams are defined in uniform rectilinear grids and hence regeneration in areas 
with little heterogeneity may be potentially wasteful. In the beam-splitting method, 
beams are automatically optimized in accordance with local heterogeneity. In other 
words, the field will be covered by minimum number of beams in a density-modulated 
manner as a result of individual independent self-similar splitting. Relative errors in 
similarity, {Nm — N)/N [m e {2,3,4}), are maximum at a; = with values —0.9%, 
— 1.3%, and —3.4%. The resultant dose errors would be smaller, due to contributions 
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of other beams, and may be tolerable. 

EfFectiveness of the splitting method was demonstrated in the experiment. The 
most prominent detour effect was the loss of range-compensated main-peak component 
in figure injc), which amounted to about 10% in dose and approximately as large as 
the distortion due to lateral particle disequilibrium. The splitting calculation and 
the measurement generally agreed well, considering that the experimental errors in 
device alignment could have been 1 mm or more. The potential sub peak for particles 
detouring around both PMMA plates was not detected, which may be natural because 
detouring itself requires scattering. The dose resolution of the MCIC system of about 
1% of the maximum should have also limited the detectability. 

In the applications to the heterogeneous phantom model, although we don't have 
reference data to compare the results with, it is natural that the splitting calculation 
with finer beams resulted in finer structures in the dose distributions. Computational 
time is always a concern in practice. In our example, the slowing factor for beam 
splitting with respect to non-splitting calculation was almost common to carbon ions 
and protons and the speed performance, a minute for 2-L volume in 1-mm grids, may 
be already acceptable for clinical applications. 

In principle, the total path length determines the computational amount for 
path integrals (jA.5p - (|A.8P and the total effective volume determines that for dose 
convolution Their Influences on the actual computational time will depend on 
algorithmic implementations (Kanematsu et al 2008a). In fact, the slowing factor for 
splitting was less than 5, which is even better than either estimation. In addition 
to common overhead that should have superficially reduced the factor, our code 
optimization with algorithmic techniques, which will be reported elsewhere, could 
have contributed to the performance. Accuracy and speed also depend strongly 
on the cutoff parameters and logical conditions in the implemented algorithm, size 
and heterogeneity of a patient model, and resolution clinically needed for a dose 
distribution. 

The automatic multiplication of tracking elements resembles a shower process 
in physical particle interactions usually calculated in MC simulations. In fact, MC 
simulations for dose calculation share many things in common. Transport and stacking 
of the elements are essentially the same and the probability for scattering may be 
equivalent to the distribution in the Gaussian approximation. As far as efficiency is 
concerned, the essential differences from the MC method are that the PB method deals 
with much less number of elements and that it does not rely on stochastic behavior of 
random numbers. 

The beam-splitting method is based on a simple principle of self-similarity and 
can be applied to any Gaussian beam model of any particle type to fill the gap 
between Monte Carlo particle simulations and conventional beam calculations in terms 
of accuracy and efficiency. However, it is difficult for beam splitting or any beam model 
in general to deal with interactions that deteriorate particle uniformity, such as nuclear 
fragmentation processes (Matsufuji et al 2005). 

5. Conclusions 

In this work, we applied our finding of self-similar nature of Gaussian distributions 
to dose calculation of heavy charged particle radiotherapy. The self-similarity enables 
dynamic, individual, and independent splitting of Gaussian beams that have been 
grown larger than the limit from lateral heterogeneity of the medium. As a result, 
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pencil beams will be arranged with optimum and modulated areal density to minimize 
overreaching and to address detouring with deflecting daughter beams. 

In comparison with a conventional calculation and a measurement, the splitting 
calculation was prominently effective in the target region with steep range adjustment 
by an upstream range compensator. The detour effect was about 10% for the 
maximum and of the same order of magnitude with lateral particle disequilibrium 
effect. In comparison between carbon ions and protons, the effects of splitting were 
not significantly different because other scattering effects were also larger for protons. 

Although performances depend strongly on physical beam conditions, clinical 
requirement, and algorithmic implementation, a typical slowing factor of the order of 
10 may be reasonably achievable for involvement of beam splitting. In fact, factor of 
5 has been achieved in our example. The principle and formulation for beam splitting 
are general and thus the feature may be added to various implementations of the PB 
algorithm in a straightforward manner. 



Appendix A. Pencil beam generation and transport 



On generation of pencil beam 6 on a plane at height zq as shown in figure [21 beam 
position rft, residual range i?^, and variances 9^1,, Oth, and t'^h are initialized as 



n,(0) = fob = fib + —vb, 

Vbz 

2 



t\{Q) 



1 

2 



ZX - Zq 



RbiO)=Ro, 

2 

cry 



Zy 



Zq Zy - Zq Sf^ 

12' 



- ^0, 
0ib{O) 



th 



-^Z^ - ZoV^Y - Zq 



(A.l) 
(A.2) 
(A.3) 



Zx Zy 

where fob = {xob,yQb, zq) is the beam-6 origin, fib — {xib,yib,0) is the beam position 
on the isocenter plane, ctx and ay are the source sizes at virtual source heights zx 
and Zy, Rq is the initial residual range, and Vb — {vbx,Vby,Vbz) is the beam direction 
radiating from the virtual sources with 



Vbz 



4+1 

4 



(A.4) 



Vbx _ _xib Vby_ _ yib 

Vbz Zx ' Vbz Zy ' 

Because nuclear interactions are effectively handled in tissue-phantom ratio D^Q{d) in 
dose calculation, number of particles Ub is modeled as invariant. 

The Fermi-Eyges theory (Eyges 1948, Kanematsu 2008c, 2009) gives increments 
of the PB parameters in step As within a density voxel by 

Afb = Vb As, ARb = ~ps As, 



0.92 



Ath^ 




Ae^b = d.oo X 10 



Aetb 



(A.5) 
(A.6) 

(A.7) 
(A.8) 



where pg and Xq/Xq^ are the effective density (Kanematsu et al 2003) and 
radiation length of the medium in units of those of water and z and m/mp are 



Dynamic splitting of pencil beams for heterogeneity corrections 



13 



the particle charge and mass in units of those of a proton. For the last physical 

step with i?t ^ and diverging AO'^b, the growth is directly given by At'^i, = 
0.0224^z-°-^^{m/mp)-°-^'^{Rb/ and then disabled by A'^b = in the unphysical 
Rb <0 region. 
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